function []=probability_variation_study(input)
%This code performs n1=20 perturbations to the tumor control probability values associated 
%to the treatment schedules included in a clinical database and performs fits to dose
%response data in terms of tumor control probability, TCP. 
%A total of n2=10 repeated optimizations are performed per probability perturbation. 
%The Logitistic TCP formulation is used with different functional dependencies 
%in the expression of the surviving fraction (SF).


%Input: One single string with input parameter values separated by spaces:
%1) tumor accelerated repopulation kickoff time in days, 2) TCP gamma parameter, 
% defining dose response curve steepness and 3) treatment site (brain or lung_).

%Input examples: '28 0.83 lung_', '7 0.7 brain'


%  A data.mat file is loaded with two clinical data arrays corresponding to a
%  NSCLC and a brain metastases RT treatment cohorts database. Two  additional
%  vectors (dep_a and dep_b) specify the functional dependencies of the 26 
%  radiobiological SF models used.


% Output: Matlab structure with:
 
%  1.	results_array: Results summary array with n1 x 26 bestfitting AICc values plus 1x26 model's
%       AICc relative standard deviation, 1x 26 model's minimum AICc value from n repeated
%       simulations, 1x26 wi values associated to the best fits from n repeated simulations,
%       1x26 mean wi values calculated for each model, 1x26 wi standard deviation values for 
%       each model, 1x26 EVRs  associated to best fits from n repeated simulations, 1x26 mean
%       EVR values calculated for each model from n repeated simulations and 1x26 EVR standard 
%       deviation for each model from n repeated simulations.

%   2.  n1*n2 Model.sim_n Matlab structures with 8 fields each: 

%       2.1	bp: 26 model' parameter values: 
%           - bp=[alphabeta a b D50 gamma lambda kickoff_time TCPmax]
%             Optimized parameters: D50, lambda and a and b (when dep_a and dep_b are not zero, respectively)
%       2.2	Loglike: 26 model's best fitting loglikelihood values. 
%       2.3	aic: 26 model's best fitting AICc values
%       2.4	prob: fitting probability associated to each point of the database for the 26 models  
%       2.5	TCPmodel: Best fitting TCP model value associated to each point of the database (26 models)
%       2.6	EQD2: EQD2 associated to the each point of the database for the 26 models
%       2.7	Effect: Treatment effect, Log(SF),associated to the each point of the database for the 26 models.
%       2.8	Evol_loglike: 26 vectors of increasingly bestfitting loglike values calculated during the fitting optimization process.

%       Structures Logit.sim_1 to Logit.sim_n2 correspond to the n2 optimizations run with the 1st dose perturbation, 
%       Structures Logit.sim_n2+1 to Logit.sim_n2+n2 correspond to the n2 optimizations run with the 2nd dose perturbation...
%       ...
%       Structures Logit.sim_(n1-1)*n2 to Logit.sim_n2*n1 correspond to the n2 optimizations run with n1th dose perturbation...

tic;

%Input parameters definition
inputs=strsplit(input);
tk=str2double(inputs(1));
gamma=str2double(inputs(2));
treatment_type=inputs(3);
%Database load
load('data.mat')
if treatment_type{1}=='brain'
    data_TCPexp=data_TCPexp_brain;
elseif treatment_type{1}=='lung_'
    data_TCPexp=data_TCPexp_lung;
end
%Output file name definition
name=strcat('Logit_', treatment_type, '_g_', num2str(gamma), '_tk_', num2str(tk),'_varp', '.mat');
% Initial parameter values: par_long0=[alpha/beta a b D50 gamma lambda_prolif t_kickoff TCPmax] 
%Simulations are run with alpha/beta=10 Gy. 
par_long0=[10 -1e-5 1e-5 25 gamma log(2)/6 tk 0.95];
%ind_opt0 = indices of the par_long0 vector corresponding to parameters free to change during the optimization process.
%Optimized parameters are [a b D50 lambda_prolif]
ind_opt0=[2 3 4 6];
T0=1;
dT=0.975;
logistic=1;
n_loop=100;
Logit.p_vector=zeros(size(data_TCPexp,1),1);
counter=1;
%Number of perturbations applied to the probability data 
rep1=20; %n1
%Number of optimizations performed per dose perturbation
rep2=10; %n2
for i=1:rep1
    data_TCPexpi=data_TCPexp;
%Tumor control probability perturbation. The number of controlled tumors is 
%sampled with binomial distribution using the number of patients/metastases 
%in the cohort and the tumor control probability, columns 1 and 3 in the database). 
%Agresti Coull estimate of tumor control probability (pAC) is recalculated using this data
%and used as perturbed probability values. 
    x2=binornd(data_TCPexp(:,1),data_TCPexp(:,3));
    pAC2=(x2+1.96^2/2)./(data_TCPexp(:,1)+1.96^2);
    data_TCPexpi(:,3)=pAC2;
    Logit.p_vector(:,end+1)=data_TCPexpi(:,3);
%The optimization process is repeated rep2 times per database dose perturbation
        for lala=1:rep2 %10
            lala
             sim=['sim' num2str(counter)];
             [Logit.(sim).bp,Logit.(sim).Loglike,Logit.(sim).aic,Logit.(sim).prob,Logit.(sim).TCPmodel,Logit.(sim).EQD2,Logit.(sim).Effect,Logit.(sim).evol_loglike]=MaxLikelihood_TCP(data_TCPexpi,par_long0,ind_opt0, T0, dT,n_loop,logistic,dep_a,dep_b,treatment_type);
             counter=counter+1;
        end
end
%Array with the 20 perturbed tumor control probability values used to fit the data
Logit.p_vector(:,1)=[];
%AIC array with 20x26 AIC values corresponding to the bestfitting models
%amongst the 10 optimizations performed for each of the 20 probability perturbations
Logit.aic=min([Logit.sim1.aic; Logit.sim2.aic; Logit.sim3.aic; Logit.sim4.aic; Logit.sim5.aic; Logit.sim6.aic; Logit.sim7.aic; Logit.sim8.aic; Logit.sim9.aic; Logit.sim10.aic]);
      for i=1:19
          Logit.aic(i+1,:)=min([Logit.(['sim' num2str(i) '1']).aic; Logit.(['sim' num2str(i) '2']).aic; Logit.(['sim' num2str(i) '3']).aic; Logit.(['sim' num2str(i) '4']).aic; Logit.(['sim' num2str(i) '5']).aic; Logit.(['sim' num2str(i) '6']).aic; Logit.(['sim' num2str(i) '7']).aic; Logit.(['sim' num2str(i) '8']).aic; Logit.(['sim' num2str(i) '9']).aic; Logit.(['sim' num2str(i+1) '0']).aic]);
      end
%Calculation of 26 Akaike weights from bestfitting models amongst the 200 optimizations performed
[Logit.wimin]=weight_calculator(min(Logit.aic));
%Calculation of 20x26 Akaike weights from the bestfitting models amongst the 10 optimizations performed per perturbation[Logit.wi]=weight_calculation(Logit.aic);
[Logit.wi]=weight_calculator(Logit.aic);

%Calculation of evidence ratios:
    for i=1:20 
        Logit.ev(i,:)=Logit.wi(i,:)./Logit.wi(i,1); 
    end
%Definition of results array:
Logit.results_array=[Logit.aic; std(Logit.aic)./mean(Logit.aic); min(Logit.aic); Logit.wimin; mean(Logit.wi); std(Logit.wi); Logit.wimin./Logit.wimin(1); mean(Logit.ev); std(Logit.ev) ];
Logit=rmfield(Logit, 'aic'); Logit=rmfield(Logit, 'wimin'); Logit=rmfield(Logit, 'ev'); Logit=rmfield(Logit, 'wi');
time=toc;
%Output file saved
save(name{1},'Logit','time');
end
